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Abstract 

Breast  cancer  is  second  only  to  lung  cancer  as  the  most  prevalent  form  of  can¬ 
cer  to  afflict  women — remaining  the  leading  cause  of  cancer  death  in  women  between 
the  ages  of  40  and  55.  Mammography  is  a  potent  weapon  in  the  fight  against  this 
lethal  disease,  due  in  large  part  to  its  widespread  availability  and  low  cost.  Despite 
the  fact  that  mammography  can  detect  small  lesions  as  early  as  two  years  before  they 
become  palpable  on  physical  exam,  between  10  and  30  percent  of  cancerous  lesions 
go  undetected  during  evaluation  by  the  radiologist.  One  approach  to  improving  de¬ 
tection  rates  involves  comparing  mammograms  of  the  same  breast  from  successive 
years.  Since  most  forms  of  breast  cancer  develop  slowly,  multiple-view  techniques 
might  be  able  to  detect  subtle  changes  indicative  of  cancerous  growth.  This  thesis 
proposes  a  computer-aided  system  designed  to  bring  two  images  into  correspondence, 
or  alignment,  so  that  they  can  be  compared  and  evaluated  for  possible  abnormalities. 
The  system  estimates  a  mapping  between  two  images  by  calculating  the  optical  flow, 
or  apparent  intensity  change,  between  a  source  and  target  mammogram.  The  effi¬ 
ciency  of  the  proposed  registration  system  is  enhanced  by  utilizing  a  multiresolution 
approach  whereby  images  are  compared  at  more  than  one  scale.  In  contrast  to  other 
registration  attempts  which  match  sets  of  morphological  features,  this  system  does 
not  require  the  identification  of  any  control  points  at  all.  This  advantage  permits 
the  system  to  perform  well  even  when  the  two  images  differ  significantly  from  one 
another.  Preliminary  results  suggest  the  potential  usefulness  of  this  system  as  part 
of  a  clinical  computer-aided  detection  (CADx)  system. 
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A  MAMMOGRAPHIC  REGISTRATION  METHOD  BASED  ON 
OPTICAL  FLOW  AND  MULTIRESOLUTION  COMPUTING 


I.  Introduction 

1.1  Introduction 

1.1.1  Background.  According  to  the  American  Cancer  Society,  breast  can¬ 
cer  is  second  only  to  lung  cancer  as  the  most  prevalent  form  of  cancer  to  afflict 
women,  but  remains  the  leading  cause  of  cancer  death  in  women  between  the  ages 
of  40  and  55.  To  put  this  problem  in  perspective,  this  year  approximately  180,000 
women  will  be  diagnosed  with  breast  cancer,  while  in  the  same  period  some  44,000 
more  women  will  lose  the  fight  against  this  deadly  disease  (60).  Computer-aided 
detection  (CADx)  technologies  hold  great  promise  in  the  effort  to  improve  early 
detection  of  breast  cancer.  This  thesis  explores  the  potential  of  analyzing  mammo¬ 
grams  taken  at  different  times  to  identify  changes  which  may  be  the  harbingers  of 
cancerous  growth. 

1.1.2  Problem  Statement.  Radiologists  typically  compare  the  current  mam¬ 
mogram  with  x-ray  images  taken  previously  in  hopes  of  detecting  changes  indicating 
possible  cancerous  tissue  growths.  Figure  1  shows  a  sequence  of  two  mammograms 
taken  at  screenings  two  years  apart.  Today  practitioners  rely,  for  the  most  part,  on 
visual  inspection  to  detect  possible  abnormal  growth  patterns.  Locating  a  region 
on  one  image  corresponding  to  a  region  of  interest  on  another  image  is  essentially  a 
judgment  call.  This  subjective  practice  is  both  time  consuming  and  error  prone. 

A  system  able  to  perform  autonomous  image-to-image  matching  would  free 
the  radiologist  to  concentrate  on  the  most  important  tasks — image  analysis  and 
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Figure  1  Time  sequence  of  mammograms,  (a)  Previous  screening  (b)  current 
screening. 


interpretation.  It  could  suggest  the  best  match  and  let  the  radiologist  decide  the 
acceptability  of  the  match.  The  ultimate  goal,  of  course,  is  to  make  the  system 
smart  enough  to  participate  in  the  change  detection  process  itself,  improving  the 
radiologist’s  ability  to  catch  the  cancer  at  a  very  early  stage. 


1.1.3  Scope.  This  thesis  will  develop  an  accurate,  reliable  method  for 
registering  a  time  series  of  mammograms.  Although  in  principle  any  number  of 
images  can  be  brought  into  a  common  alignment,  this  thesis  will  concentrate  on 
a  sequence  of  two  images.  Specifically,  the  matching  will  be  constrained  to  two 
images  of  the  same  breast  taken  at  different  times.  The  registration  method  will  not 
support  bi-lateral  comparisons  of  the  left  and  right  breasts.  The  algorithm  will  be 
based  solely  on  optical  flow  approaches. 

The  registration  scheme  will  interface  seamlessly  with  other  image  analysis 
activities.  The  registration  procedure  is  not  an  end  in  itself;  rather,  it  is  merely 
one  of  many  preprocessing  activities.  The  driving  force  behind  the  development  of 
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registration  techniques  is  the  need  to  detect  cancer-induced  tissue  changes  as  early 
as  possible.  Comparing  mammograms  from  the  current  and  previous  screenings  has 
the  potential  to  detect  warning  signs  long  before  they  would  be  discovered  by  single¬ 
image  analysis.  With  this  goal  in  mind,  the  registration  technique  should  enhance 
post-registration  mammogram  analysis  as  much  as  possible. 

1.1.4  Methodology.  Although  a  series  of  mammograms  taken  from  the 
same  view  (i.e.  craniocaudal  or  mediolateral)  exhibit  many  similarities  over  time, 
the  images  may  also  appear  quite  different.  Certain  types  of  distortion  may  be  traced 
to  the  imaging  equipment,  while  others  can  be  attributed  to  the  image-taking  pro¬ 
cedure  itself.  An  example  of  equipment-based  distortion  is  variation  in  background 
illumination,  causing  one  image  to  appear  brighter  than  the  other.  Differences  in 
breast  compression  from  one  session  to  the  next  cause  its  overall  shape  to  appear 
different.  Although  care  is  taken  during  the  mammography  session  to  ensure  proper 
alignment  of  the  breast,  there  is  very  little  the  technician  can  do  to  ensure  exact 
positioning  of  the  breast. 

The  complexity  of  the  task  makes  it  clear  that  there  is  no  silver  bullet  solution. 
Instead,  the  synergistic  effect  of  several  techniques  may  make  the  solution  more 
tractable.  The  method  developed  here  is  based  on  a  series  of  processing  stages,  each 
trying  to  further  constrain  the  problem.  In  particular,  the  method  follows  a  coarse- 
to-fine  approach  in  which  global  matching  precedes — and  indeed  enables — matching 
at  finer  and  finer  levels  of  resolution.  This  approach  mitigates  the  potential  for  the 
solution  to  be  less  than  acceptable. 

1. 1.4.I  Global  Alignment.  The  first  problem  encountered  has  to  do 
with  the  orientation  problem  discussed  above.  Clearly,  the  two  breasts  in  the  mam¬ 
mograms  must  first  be  aligned  to  the  same  relative  position  and  orientation  before 
any  serious  attempt  can  be  made  to  correspond  regions  within  the  breast  struc¬ 
ture.  A  simple  and  effective  way  to  ensure  optimal  alignment  involves  identifying 
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an  appropriate  affine  transformation  from  one  image  to  the  other.  An  affine  trans¬ 
formation  can  be  considered  a  parametric  estimate  of  global,  linear  differences  with 
respect  to  rotation,  translation,  and  scale  (25).  The  method  proposed  here  solves 
for  an  optimal  set  of  parameters  using  a  genetic  algorithm  and  realigns  the  source 
image  to  the  target  image  based  on  these  values. 

1.1. 4. 2  Regional  Deformation.  The  affine  transformation  carried  out 
in  the  previous  step  does  not  correct  nonlinear  deformations  such  as  compression 
and  dilation.  The  approach  taken  here  involves  identifying  the  apparent  motion 
(directional  change  in  intensity)  between  the  source  and  target  images.  The  proposed 
method  uses  the  optical  flow  as  a  measure  of  change  between  images.  The  inter-image 
velocity  information  is  then  used  to  deform  the  source  image  into  an  intermediate 
form  more  closely  resembling  the  target  image.  The  transformation  is  carried  out  by 
an  appropriate  interpolation  function. 

1.1. 4- 3  Local  Similarity  Assessment.  During  the  registration  pro¬ 
cess,  a  similarity  measure  between  the  transformed  source  image  and  the  target 
image  must  be  evaluated.  In  this  system,  a  neural  network  approach  based  on 
the  Kohonen  Self-Organizing  Map  (SOM)  is  taken  to  develop  a  problem-dependent 
similarity  metric.  The  network  learns  to  segment  an  image  by  partitioning  it  into 
homogeneous  regions  based  on  the  spatial  distribution  of  image  intensity  values.  A 
traditional  subtraction  method  is  also  used  to  measure  the  similarity  between  two 
images. 

1.1.5  Overview.  The  remaining  chapters  are  organized  as  follows.  Chapter 

II  discusses  previous  and  on-going  work  in  the  area  of  image  registration.  Chapter 

III  describes  the  proposed  mammogram  registration  system.  Experimental  results 
and  conclusions  are  discussed  in  Chapters  IV  and  V,  respectively. 
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II.  Background 

This  chapter  first  motivates  the  need  for  registering  mammograms  from  a  health  care 
perspective.  It  then  discusses  several  technical  aspects  of  image  processing  relevant 
to  the  registration  problem.  Finally,  the  major  issues  involved  in  registration  are 
discussed,  highlighting  key  research  in  each  area. 

2.1  Breast  Cancer 

Breast  cancer  is  second  only  to  lung  cancer  as  the  leading  form  of  cancer  to 
strike  women  (63).  While  mammography  has  proven  to  be  a  highly  effective  tool  in 
the  early  detection  of  breast  cancer,  between  10  and  30  percent  of  cases  diagnosed 
as  cancer  had  negative  mammogram  interpretations. 

2.2  Computer-Aided  Diagnosis  (CADx) 

Considerable  effort  is  being  devoted  to  the  development  of  automated  systems 
capable  of  assisting  the  radiologist  in  mammogram  analysis  and  interpretation  (50, 
37,  61).  A  typical  CADx  system  uses  sophisticated  image  processing  techniques  to 
detect  either  microcalcifications  or  masses  which  might  otherwise  go  unnoticed.  It 
is  important  to  note  that  these  computer-based  systems  are  designed  to  assist — not 
replace — the  trained  radiologist. 

Due  to  the  complexity  of  mammogram  interpretation,  radiologists  base  their 
decisions  on  all  available  information.  In  practice,  this  means  comparing  all  four 
views  (craniocaudal  and  mediolateral  oblique  for  each  breast)  from  the  current 
screening  and  prior  screenings.  Since  most  forms  of  breast  cancer  develop  slowly, 
comparison  of  mammograms  taken  at  different  times  can  provide  evidence  of  a 
growth  not  previously  detected  (40,  77). 

Although  most  CADx  implementations  currently  under  development  or  in  clin¬ 
ical  test  are  limited  to  single  image  analysis,  attempts  are  underway  to  develop 
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multi-view  capable  systems  (50).  Dinten  et  al.  (18)  propose  a  system  which  com¬ 
bines  analysis  from  single  mammograms  with  a  bilateral  comparison  between  left 
and  right  mammograms.  The  primary  stumbling  stock  to  the  development  of  these 
CADx  systems  is  the  need  to  register  two  images  before  they  can  be  compared  and 
analyzed.  The  variation  between  two  same-breast  mammograms  taken  at  different 
times  can  be  enormous,  rendering  the  registration  task  impractical.  For  example, 
differences  in  breast  positioning  and  compression  can  result  in  appreciable  distortion. 
Until  these  complications  have  been  resolved,  progress  in  multi-view  CADx  system 
development  will  remain  slow. 

2.3  Image  Processing  Methods 

Image  processing  is  often  divided  into  three  stages:  preprocessing,  feature  ex¬ 
traction,  and  classification  (59).  Preprocessing  consists  of  low-level  techniques  de¬ 
signed  to  transform  the  input  data  into  a  form  useful  to  the  other  stages.  Feature 
extraction  identifies  and  isolates  characteristics  (patterns)  in  the  preprocessed  data 
based  on  semantic  merit.  Finally,  extracted  features  are  combined  and  manipulated 
so  that  the  original  problem  can  be  solved.  The  following  subsections  examine  some 
representative  image  processing  methods  most  important  to  solving  the  registration 
problem. 

2.3.1  Noise  Reduction.  The  useful  information  conveyed  by  images  is  par¬ 
tially  degraded  by  unwanted  signals  in  the  form  of  noise.  Origins  of  noise  include  the 
image  capturing  device  (e.g.  mammography  unit),  the  subsequent  film  development 
process,  and  the  digitization  of  this  film  for  computer  analysis.  The  presence  of  noise 
may  corrupt  the  image  and  lessen  the  radiologist’s  ability  to  detect  an  abnormality 
in  the  mammogram  (39). 

Noise  reduction  is  most  often  achieved  by  smoothing  the  image  with  one  or 
more  filters  (28).  Although  Gaussian  filters  have  found  wide  support  in  the  image 
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analysis  community,  they  have  the  undesirable  property  of  smoothing  over  discontin¬ 
uous  as  well  as  continuous  areas,  possibly  removing  valuable  boundary  information 
in  the  process  (2).  Saint-Marc  et  al.  propose  an  adaptive  smoothing  technique  de¬ 
signed  to  retain  the  benefits  of  Gaussian  smoothing  while  preserving  discontinuities 
(60). 


2.3.2  Image  Segmentation.  Segmentation  refers  to  the  process  of  parti¬ 
tioning  an  image  into  homogeneous  regions  based  on  texture,  intensity,  or  any  other 
problem-dependent  measure.  Techniques  for  segmenting  an  image  can  be  grouped 
into  three  general  classes:  edge  detection,  thresholding,  and  region  growing  (25). 
An  important  application  of  this  technique  involves  extracting  the  breast  structure 
from  the  background  image  (5,  13,  34,  48,  53,  69). 

2.4  Registration  Methods 

Registration  is  the  process  of  establishing  the  best  possible  correspondence 
between  two  objects — images  in  this  case.  Brown  provides  an  exhaustive  survey 
of  registration  techniques  and  offers  an  excellent  taxonomy  for  making  sense  of  the 
many  approaches  available  (9).  The  following  discussion  on  registration  methods 
and  characteristics  will  follow  Brown’s  organizational  framework. 

2.4.I  Correlation.  Correlation-based,  sometimes  called  block  matching, 
techniques  attempt  to  register  two  images  using  similarity  measures.  As  Brown 
points  out,  correlation  matching  is  more  a  similarity  metric  than  a  registration 
method  unto  its  own.  Combined  with  template  matching,  cross-correlation  can  pro¬ 
vide  a  useful  measure  of  how  well  a  region  in  the  source  image  spatially  matches  a 
region  in  the  target  image.  In  the  standard  block  matching  method  given  by 

«({<*.,  »<})  =  E  IIAMi)  +  (1) 

peR 
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the  error  E({ui,Vi})  is  minimized  over  a  small  image  region  where 
represents  the  optimum  displacement  vector  (25). 

Correlation  is  particularly  useful  when  spatial  structure  remains  relatively  con¬ 
stant  between  two  images.  For  example,  two  different  aerial  photographs  taken  of 
the  same  land  area  fit  this  criterion.  A  translational  shift  in  one  image  relative  to  the 
other  will  have  no  effect  on  the  structure  of  a  given  region.  However,  a  relative  dif¬ 
ference  in  scale  between  the  two  images  will  introduce  considerable  noise  (error)  into 
the  correlation  measurement.  Caves  (12)  demonstrates  the  effectiveness  of  template 
matching  in  registering  two  land  mass  images  taken  by  Synthetic  Aperture  Radar 
(SAR). 

2.4-2  Fourier  Methods.  Whereas  correlation  is  based  on  spatial  charac¬ 
teristics,  Fourier  methods  operate  in  the  frequency  domain.  Calway  (11)  presents 
a  multiresolution  Fourier  Analysis  approach  to  the  problem  of  estimating  dispar¬ 
ity  and  motion  in  image  registration.  This  method  is  appropriate  if  the  frequency 
components  of  the  image  correspond  to  the  actual  disparity. 

2.4.3  Point  Mapping.  Mapping  a  set  of  points  in  a  source  image  to  its 
counterpart  in  a  target  image  is  a  very  powerful  registration  method  if  stable  corre¬ 
spondence  points  can  be  extracted  from  each  image  (54).  However,  significant  vari¬ 
ations  between  the  two  images  can  make  these  fiducial  markers,  or  control  points, 
very  difficult  to  identify. 

This  approach  has  recently  been  applied  to  the  registration  of  two  mammo¬ 
grams  with  considerable  success  (72).  It  should  be  noted  that  the  images  presented 
in  (72)  have  easily  identifiable  fiducial  markers — an  exception  rather  than  the  rule  for 
most  mammogram  comparisons.  Bookstein  (8)  has  shown  how  this  point-to-point 
mapping  can  be  decomposed  into  its  constituent  linear  and  nonlinear  components. 
This  type  of  analysis  is  useful  in  understanding  the  displacement  forces  responsible 
for  transforming  one  image  into  the  other. 


8 


2-4-4  Deformable  Models.  Deformable  models  are  attracting  considerable 
attention  in  many  areas  of  image  processing  and  computer  analysis.  Durbin  (20) 
demonstrated  the  power  of  elastic  nets  by  finding  a  good  solution  to  the  NP-complete 
Traveling  Salesman  Problem.  Also  known  as  snakes,  elastic  nets,  and  active  contour 
models,  these  versatile  tools  offer  a  fresh  perspective  to  long-standing  challenges  such 
as  scene  segmentation,  boundary  extraction,  and  image  registration.  Many  of  the 
approaches  under  investigation  can  trace  their  roots  to  the  snake  model  developed  by 
Kass  et  al.  (36).  These  investigators  used  this  one-dimensional  deformable  model 
to  seek  out  curves  in  images.  The  snake  is  really  a  spline  which  reaches  its  desired 
configuration  by  minimizing  a  three-component  energy  function: 


Dsnake  —  J  Esna}ie {vi^sfjds  —  j  -S|nt(v(s))  “t"  Eimage (v (s) )  +  Econ (v (3) )c?s  (2) 

The  first  component,  the  internal  energy  term,  represents  the  energy  stored 
in  the  deformable  model  due  to  bending.  The  second  term,  the  external  energy, 
corresponds  to  image  features.  The  third  component  of  the  energy  function  allow 
the  user  to  add  additional  constraints  to  the  solution. 

A  number  of  investigators  have  applied  deformable  models  to  difficult  problems 
in  medical  image  analysis.  Lobregt  et  al.  (45)  use  an  active  contour  model  to  locate 
arteries  in  vascular  x-ray  images  and  tumors  in  MR  images.  Ranganath  (55),  Cohen 
et  al  (14),  and  Coppini  et  al.  (15)  apply  similar  approaches  to  the  automatic 
extraction  of  left  ventricle  boundaries  from  cardiac  imagery.  Yue  et  al.  (75)  detect 
rib  borders  with  a  snake.  A  deformable  surface  model  is  used  by  Sandor  et  al.  to 
match  MR  brain  images  with  a  brain  atlas.  Davatzikos  and  Prince  (16)  developed 
a  technique  for  mapping  the  outer  cortex  in  brain  images  using  deformable  ribbons. 


9 


2.5  Registration  Considerations 

2.5.1  Feature  Selection.  The  image  features  selected  to  guide  the  regis¬ 
tration  process  can  dramatically  affect  the  final  outcome.  Simply  stated,  features 
which  represent  the  image  well  should  produce  a  more  accurate  registration  result 
than  features  less  representative  of  the  images.  Features  can  be  categorized  as  either 
structural  and  statistical  in  nature.  Structural  features  represent  morphological  as¬ 
pects  of  the  image,  such  as  edges  and  landmarks.  Statistical  features,  on  the  other 
hand,  capture  textural  properties  of  the  image  in  the  form  of  intensity  variance,  mo¬ 
ments,  etc.  Structural  features  are  sometimes  called  intrinsic  features  because  they 
are  derived  directly  from  the  image. 

2.5.2  Similarity  Measure.  Closely  related  to  feature  selection  is  the  decision 
of  how  to  measure  the  similarity  between  two  images.  Anandan  (1)  uses  a  simple 
sum  of  squared  difference  (SSD)  measure,  given  by  the  equation 


E({ui,  Vi})  =  52  [h(xi  +  uh yi  +  Vi)  -  I0(xi,  yi)]2  (3) 

t 

to  measure  the  difference  between  two  image  frames.  This  optimization  equa¬ 
tion  can  be  readily  solved  by  a  least  squares  or  other  iterative  numerical  method. 

Rogowska  (56)  uses  the  normalized  cross-correlation  equation 


NCORn  =  <*«)  (4) 

y  [C^*i  —  HAij)  •  (A-ij  —  fJ^Aij)]  [(-^  —  Hr)  •  (U  —  Hr)} 

to  perform  a  type  of  dynamic  image  analysis  called  similarity  mapping  between 
two  consecutive  frames  of  a  dynamic  MRI  temporal  sequence.  The  author  presents 
evidence  that  this  technique  can  identify  different  image  structures  having  clinical 
interest. 
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2.5.3  Search  Models.  Having  selected  a  set  of  features  and  similarity  mea¬ 
sure,  a  search  strategy  appropriate  to  the  type  of  registration  problem  must  be 
identified.  Specifically,  the  a  transformation  model  appropriate  to  the  registration 
problem  at  hand  must  be  selected.  The  model  should  match  the  expected  transfor¬ 
mation  as  closely  as  possible.  Transformation  models  typically  used  by  registration 
procedures  are  fall  into  one  of  four  families:  rigid,  affine,  projective,  or  curved  (70). 

2.5.3. 1  Affine  Transformations.  The  simplest  types  of  transforma¬ 
tions  are  affine.  An  affine  transformation  is  restricted  to  changes  resulting  from 
translation,  rotation,  and  shear.  Therefore,  lines  parallel  before  the  transformation 
will  remain  parallel  in  the  transformed  image.  The  parametric  affine  model  defined 
mathematically  by 


u(x,  y‘,  a) 


sxcosOx  —SySindy 

X 

+ 

dx 

sxsindx  Sycosdy 

. v 

dy 

(5) 


will  transform  each  point  in  the  source  image  to  its  corresponding  location  in 
the  target  image.  Fuh  (25)  use  an  affine  model  to  estimate  the  velocity  field  between 
two  frames  of  an  aerial  image  sequence. 


2.5.4  Search  Strategy.  The  final  component  to  make  the  registration  sys¬ 
tem  complete  is  the  search  strategy.  This  decision  is  critical  because  it  will  decide 
important  issues  such  as  registration  accuracy  and  computational  resource  require¬ 
ments. 


2.5.4-1  Algorithm-based.  Many  traditional  search  strategies  are  based 
on  well-defined  algorithms.  If  the  registration  problem  can  be  subdivided  into  smaller 
and  smaller  subproblems,  the  matching  problem  may  be  amenable  to  a  dynamic 
programming  approach.  However,  the  exhaustive  nature  of  dynamic  programming 
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algorithms  can  cause  computation  time  and  resource  requirements  to  become  pro¬ 
hibitive. 


2. 5.4 - 2  Relaxation-based.  Several  researchers  have  developed  opti¬ 
mization  approaches  to  solve  the  registration  problem.  Relaxation  methods  are  typ¬ 
ically  used  to  propagate  local  information  throughout  the  system,  enabling  a  more 
accurate  global  solution.  In  general,  these  methods  search  for  an  optimal  solution  by 
minimizing  an  objective  energy  function.  Comprehensive  discussions  on  relaxation 
techniques  and  their  application  to  image  processing  can  be  found  in  (32,  57,  58). 

2.5. 4- 3  Deterministic  Relaxation.  Numerical  analysis  theory  offers 
a  number  of  techniques  to  solve  optimization  problems  through  deterministic  relax¬ 
ation  (47).  Iterative  methods,  such  as  the  Jacobi,  Gauss-Seidel,  and  simultaneous 
over-relaxation  (SOR)  variants,  gradually  converge  to  a  stable  minimum  (under  the 
appropriate  conditions).  Recently,  the  Hopfield  network  has  been  used  to  perform 
the  optimization  task  (13,  17,  32,  41,  43,  53,  54,  68,  78).  The  Hopfield  network 
is  a  fully-connected  single-layer  neural  network  which  can  minimize  a  cost  function 
representing  the  problem  at  hand.  This  architecture  lends  itself  to  parallel  computer 
implementations . 

2.5. 4.4  Stochastic  Relaxation.  Complementing  deterministic  relax¬ 
ation  techniques  are  those  based  on  stochastic  processes  (24,  26,  49,  69).  Geman  and 
Geman  (26)  demonstrate  how  image  restoration  problems  can  be  addressed  with 
a  probabilistic  relaxation  technique  called  simulated  annealing.  This  approach  can 
avoid  local  minima  by  allowing  states  to  be  accepted  that  actually  increase  the  sys¬ 
tem’s  energy.  The  probability  that  such  an  energy-increasing  state  will  be  tolerated 
is  determined  by  an  annealing  schedule.  Initially  set  to  a  high  value,  the  acceptance 
probability  (based  on  a  Gibbs  distribution)  is  gradually  lowered  until  finally  only 
energy-decreasing  states  are  allowed.  Although  very  effective  at  reaching  an  opti- 
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mal  global  minimum,  simulated  annealing  methods  generally  suffer  from  very  slow 
convergence  rates  compared  to  other  deterministic  relaxation  methods. 

2.5. 4-5  Genetic  Algorithms.  Genetic  algorithms  offer  a  biologically- 
motivated  alternative  to  traditional  search  approaches.  Patterned  after  DNA  repli¬ 
cation  in  living  creatures,  genetic  algorithms  conduct  global  searches  with  the  help 
of  mutation  and  crossover  of  chromosomes  (27,  30).  These  artificial  chromosomes 
represent  partial  solutions  to  the  larger  problem  to  be  solved.  Over  the  course  of 
several  generations,  these  building  blocks  are  combined  so  as  to  evolve  an  optimal 
solution.  Genetic  algorithms  thrive  on  search  spaces  too  messy  to  be  solved  by  more 
structured  search  strategies. 

2. 6  Optical  Flow 

The  analysis  of  image  sequences  requires  an  understanding  of  how  they  differ 
from  one  another.  In  the  context  of  image  registration,  we  wish  to  map  a  set  of 
features  in  the  source  image  to  the  corresponding  set  of  features  in  the  target  image. 
These  features  can  be  intrinsic  image  structures  like  contours,  or  they  can  be  statis¬ 
tical  measures  such  as  contrast,  moments,  etc.  Although  several  techniques  exist  to 
estimate  this  interimage  mapping,  one  approach  which  has  proven  extremely  useful 
is  based  on  the  estimation  of  optical  flow.  Luettgen  et  al.  define  optical  flow  as  the 
”  apparent  velocity  vector  field  corresponding  to  the  observed  motion  of  brightness 
patterns  in  successive  image  frames”  (46). 

Barron  et  al.  (3)  provide  a  comprehensive  survey  of  optical  flow  approaches 
and  assess  nine  of  the  most  popular  methods  in  use  today.  The  authors  place  optical 
flow  estimation  techniques  into  four  categories:  differential,  region-based  matching, 
energy-based,  and  phase-based.  Differential  methods  use  spatio-temporal  deriva¬ 
tives  to  estimate  a  dense  velocity  field  relating  the  source  image  to  the  target  image. 
Region-based  techniques,  on  the  other  hand,  rely  on  similarity  measures  between 
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features  in  the  image  sequence  to  derive  a  sparse  optical  flow  field.  The  third  cate¬ 
gory,  energy-based  optical  flow  techniques,  apply  velocity-tuned  filters  to  produce  a 
flow  estimate.  In  a  similar  vein,  phase-based  techniques  use  the  outputs  of  band-pass 
filters  to  derive  the  flow  estimate. 

2.6.1  Horn  and  Schunck  Method.  In  many  types  of  image  matching,  this 
set  of  consistent  features  is  difficult,  if  not  impossible,  to  identify.  Horn  and  Schunck 
(31)  established  the  basis  for  the  development  of  homogeneous  differential  optical 
flow  methods  based  on  image  gradient  information.  Their  method  is  predicated  on 
the  brightness  constancy  equation 

dE  ,  . 

—  =  Exu  +  Eyv  +  Et  =  0  (6) 

where 


b,  _  dE(x,y,t)  „  _dE(x,y,t )  „  _  dE(x,y,t) 
Ex=  d~x  'Ev=  d~y  ’ Et=  dt 

represent  the  image  spatial  and  temporal  intensity  derivatives  and 


(u,v)  =  (dx/dt,dy/dt)  (8) 

denotes  the  derived  optical  flow.  The  brightness  constraint  assumes  that 
”  changes  in  image  brightness  are  due  only  to  motion  in  the  image  frame”  (46). 

This  single  linear  equation  in  two  unknowns  results  in  an  ill-posed  problem, 
requiring  more  information  for  an  acceptable  solution  (4).  An  ill-posed  problem 
denotes  a  solution  that  is  in  general  neither  unique  nor  stable  (76).  One  well- 
established  technique  to  condition  this  type  of  underconstrained  problem  is  to  add  a 
regularization  term.  Poggio  et  al.  (52)  describe  how  regularization  techniques  can  be 
used  to  tame  a  variety  of  difficult  problems,  including  optical  flow  calculations.  The 
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regularization  term  is  not  constrained  to  assume  a  particular  form;  it  simply  needs 
to  add  sufficient  information  to  make  the  problem  well-posed.  Horn  and  Schunck 
(31)  attack  the  problem  by  incorporating  a  smoothness  constraint,  resulting  in  a 
minimization  problem  of  the  form 

$=  /  [  (Exu  +  Eyv  +  Et)2  +  a2(ul  +  u2  +  vl  +  vl)dxdy  (9) 

J  Jlmage 

where  a  represents  the  regularization  parameter  (smoothness,  in  this  case). 

This  quadratic  variational  functional  can  be  recast  into  a  set  of  Euler-Lagrange 
equations  of  the  form 


( Exu  +  Eyv  +  Et)Ex  =  a2  A  u(Exu  +  Evv  +  Et)Ey  =  a2Av  (10) 

2. 6.2  Alternative  Methods.  Several  researchers  have  proposed  modifica¬ 
tions  to  the  basic  Horn  and  Schunck  framework  in  an  attempt  to  improve  accu¬ 
racy  of  the  optical  flow  calculation.  Black  et  al.  (6,  7)  propose  a  robust  estima¬ 
tion  framework  designed  to  minimize  the  negative  effect  of  violations  to  the  basic 
assumptions — brightness  constancy  and  spatial  smoothness.  Specifically,  they  incor¬ 
porate  p-functions,  such  as  the  Lorentzian  and  Geman-McClure  variants,  to  atten¬ 
uate  the  influence  of  outlying  measurements  on  the  solution.  Optical  flow  is  then 
estimated  by  minimizing  the  objective  function 


E 


ADM(VJ(*Mx) +/,(*)), +  ■£  pd|«(*)-«(*)||,<r») 

zeG(x ) 


(11) 
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The  first  term  enforces  fidelity  to  image  intensity  values  as  determined  by  the 
control  parameter  an,  while  the  second  term  regulates  smoothness  as  determined  by 
another  control  parameter  as  and  serves  as  the  regularization  parameter. 

Fleet  (23)  uses  local  phase  information  to  calculate  optical  flow,  while  Fuh 
(25)  extracts  regional  optical  flow  estimates  from  the  observed  affine  transformation 
in  each  region.  Nagel  and  Enkelmann  provide  an  excellent  discussion  on  the  role  of 
smoothness  constraints  in  optical  flow  formulations  (51).  Fan  et  al.  (22)  describe 
an  image  motion  estimation  method  based  on  the  Expectation-Maximization  (EM) 
technique  and  affine  transformations.  Verri  (71)  also  gives  a  good  discussion  on 
motion  fields  and  optical  flow  computation.  Zhang  and  Hanauer  (76)  apply  the 
concepts  of  mean  field  theory  to  address  the  optical  flow  problem. 

2.7  Multiscale  Methods 

The  notion  of  multiresolution  image  analysis  is  not  new.  Research  efforts  into 
new  analysis  approaches  integrating  information  generated  on  more  than  one  scale 
have  been  reported  since  the  mid  1970’s  (58).  Burt  (10)  introduced  the  idea  of 
a  Laplacian  pyramid  to  create  a  compact  image  code.  Although  he  developed  this 
technique  to  improve  data  transmission  performance,  this  multiscale  analysis  tool 
has  found  application  in  the  image  analysis  realm  and  has  been  applied  to  hyper- 
spectral  image  data  fusion  (73).  The  Gaussian  kernel  is  generally  preferred  because 
”it  is  the  only  kernel  that  leads  to  monotonic  destruction  of  detail  under  repeated 
blurring... no  new  details  will  appear  at  a  given  scale  s  that  did  not  exist  at  smaller 
scales”  (44).  Repeated  blurring  of  the  source  image  defines  a  continuously-generated 
family  of  images  referred  to  as  scale  space. 

Multiscale  methods  offer  levels  of  analysis  and  manipulation  not  generally  avail¬ 
able  using  single-scale  methods.  Specifically,  they  employ  coarse-to-fine  strategies 
which  decompose  an  image  (or  sequence  of  images)  into  more  than  one  level  of 
resolution.  This  method  is  extremely  useful  in  image  registration  because  motion, 
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or  intensity  change  between  two  images,  will  often  be  present  to  different  degrees. 
In  other  words,  there  might  be  rapid  change  at  the  breast  edge  along  with  much 
smaller  local  change  from  region  to  region.  Simple  optical  flow  methods  (e.g.  Horn 
and  Schunck)  assume  a  linear  intensity  change  within  the  measurement  area,  but 
this  assumption  is  often  violated. 

Multiscale  methods  address  this  well-known  ’’aperture  problem”  by  looking 
at  the  problem  from  different  levels  of  resolution.  Using  this  approach,  small  in¬ 
tensity  changes  can  be  investigated  at  fine  levels  while  larger  intensity  changes  can 
be  explored  at  an  appropriate  level  of  resolution.  For  the  purpose  of  optical  flow 
calculations,  we  look  for  the  level  at  which  a  given  intensity  change  is  linear.  This 
hierarchical  decomposition  ensures  that  all  degrees  of  motion  are  analyzed  at  the 
level  appropriate  to  their  magnitude. 

Several  investigators  have  developed  optical  flow  estimation  methods  based 
on  this  multiresolution  approach  (4,  21,  33,  62,  65,  66,  64,  67,  74).  Szeliski  et  al. 
(65,  66,  64)  base  their  methods  on  a  hierarchical  family  of  splines.  Battiti  et  al. 
(4)  address  the  shortcomings  of  homogeneous  differential  methods  (e.g.  Horn  and 
Schunck)  by  estimating  the  relative  error  associated  with  each  optical  flow  calcula¬ 
tion.  An  optical  flow  measurement  with  error  exceeding  a  specified  threshold  will  be 
recalculated  at  the  next  finer  level,  while  those  below  the  threshold  (high  confidence 
estimates)  axe  retained.  In  a  similar  vein,  Enkelmann  (21)  and  Hwang  (33)  also 
estimate  the  relative  error;  their  approach  incorporates  the  error  estimate  directly 
into  the  next  calculation  iteration.  Hwang  and  Lee  (33)  apply  multigrid  techniques 
to  their  multiresolution  pyramid  with  excellent  results.  Luettgen  et  al  (46)  de¬ 
velop  a  multiscale  regularization  approach  for  robust  optical  flow  estimation.  The 
performance  of  such  multiresolution  frameworks  can  be  dramatically  improved  by 
incorporating  true  multigrid  algorithms  that  make  explicit  use  of  interlevel  residual 
errors  to  speed  up  the  relaxation  process.  Douglas  and  Terzopoulos  (19,  67)  present 
a  comprehensive  discussion  of  multigrid  techniques  for  improving  convergence  rates. 
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2.8  Summary 

Mammography  has  proven  to  be  a  powerful  tool  in  the  fight  against  breast 
cancer.  When  interpreting  mammograms,  radiologists  normally  compare  the  current 
mammogram  set  with  those  taken  previously.  Researchers  are  actively  developing 
computer-aided  detection  (CADx)  systems  capable  of  helping  the  radiologist  make 
more  informed  decisions.  Much  of  today’s  CADx  research,  however,  still  focuses 
on  the  analysis  of  one  image  at  a  time.  Although  single-view  analysis  is  extremely 
important,  a  sequence  of  mammograms  contains  information  about  the  growth  of 
cancer  over  time.  Since  most  types  of  breast  cancers  grow  slowly,  the  simultaneous 
analysis  of  mammograms  spanning  more  than  one  year  may  significantly  improve 
the  detection  of  cancers  too  difficult  to  identify  on  a  single  image. 

A  key  technology  in  multiple-image  processing  is  registration — matching  loca¬ 
tions  on  one  image  with  the  corresponding  locations  on  the  other.  The  significant 
differences  typically  encountered  between  two  images  often  make  registration  very 
difficult.  In  the  case  of  mammography,  it’s  a  non-linear  problem  requiring  a  non¬ 
linear  solution.  Using  optical  flow  techniques  to  estimate  the  disparity  between  two 
images  and  facilitate  registration  is  not  a  new  idea.  However,  these  powerful  tech¬ 
niques  have  not  as  yet  been  applied  to  the  registration  of  mammograms.  Optical 
flow  techniques  offer  a  potentially  powerful  approach  to  this  complex  problem.  This 
thesis  develops  a  registration  method  based  on  optical  flow  and  explores  both  the 
advantages  and  disadvantages  of  this  approach. 
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III.  Methodology 


This  chapter  describes  the  proposed  mammogram  registration  system. 

3.1  System  Overview 

The  mammogram  registration  system  developed  in  this  thesis  consists  of  sev¬ 
eral  modules  each  designed  to  perform  a  specific  function.  It  relies  on  both  ana¬ 
lytical  techniques,  such  as  iterative  relaxation,  where  appropriate,  and  appeals  to 
biologically-motivated  methods,  such  as  genetic  algorithms  and  neural  networks, 
where  they  outperform  the  traditional  approaches.  It  is  based  on  a  multiresolu¬ 
tion  hierarchy  to  facilitate  analysis  at  more  than  one  scale.  This  approach  enables 
the  system  to  deal  with  registration  and  tumor  detection  at  the  most  appropriate 
level(s)  of  resolution.  From  a  philosophical  standpoint,  this  system  is  designed  to 
be  an  integrated  approach  to  both  mammogram  registration  and  tumor  detection. 
Specifically,  the  proposed  system  uses  velocity  field  estimates  between  two  images  to 
simultaneously  bring  two  images  into  correspondence  and  identify  potential  tumor 
growth  areas.  Figure  2  gives  a  bird’s-eye  view  of  the  entire  registration  system.  Fig¬ 
ure  3  displays  a  time  series  of  four  mammograms  for  the  same  woman  taken  during 
different  screenings.  The  registration  system  must  align  these  images  to  help  the 
radiologist  compare  corresponding  regions  in  hopes  of  detecting  a  growth  if  one  is 
present. 

The  preprocessing  module  is  responsible  for  low-level  data  manipulation  to  en¬ 
sure  both  images  are  as  similar  as  possible.  The  global  alignment  module  reduces 
differences  in  position  and  orientation  prior  to  estimating  intensity  differences  be¬ 
tween  two  images.  The  segmentation  module  extracts  each  breast  structure  from 
the  background  to  help  define  for  the  system  what  to  register  and  what  to  ignore. 
The  multiscale  pyramid  generator  decomposes  each  image  into  a  specified  number 
of  scales.  The  number  of  levels  is  not  fixed  and  can  be  adjusted  to  meet  the  current 
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MAMMOGRAM  REGISTRATION  SYSTEM 


Preprocess 


Global  Alignment 


Segment  From 
Background 


Construct  Pyramid 


Compute  Optical 
Flow 
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Figure  2  Proposed  mammogram  registration  system. 


20 


a006b00b.mam 


a006bl0u.mam 


50  100  150  200  250  50  100  150  200  250 


(c)  (d) 

Figure  3  Time  sequence  of  mammograms  taken  during  four  consecutive  screen¬ 
ings.  (a)  First  screening  (b)  second  screening  (c)  third  screening  (d)  last 
screening. 
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need.  The  heart  of  the  registration  system  is  the  optical  flow  computation  module. 
This  module  estimates  the  apparent  difference  between  the  two  images  being  reg¬ 
istered.  The  warping  module  takes  the  output  from  the  optical  flow  computation 
module  and  generates  a  warped  image  based  on  the  apparent  motion  from  source 
to  target  image.  Finally,  the  evaluator  module  measures  the  similarity  between  the 
warped  image  and  target  image  and  decides  whether  the  sequence  must  be  repeated. 
The  similarity  maps  generate  can  be  used  to  evaluate  the  possibility  of  abnormal 
growth  in  the  second  mammogram. 

3.2  Preprocessing 

3.2. 1  Noise  Reduction.  Noise  due  to  film  granularity  and  introduced  by  the 
digitization  process  will  always  be  present  in  mammogram  images.  Median  filtering 
has  been  shown  to  be  an  effective  tool  in  removing  this  type  of  noise.  In  this  system, 
images  are  filtered  using  a  median  filter  with  a  5  x  5  kernel. 

3.2.2  Normalization.  The  first  problem  encountered  by  any  image  registra¬ 
tion  process  is  the  variety  of  differences  which  may  be  present  in  the  images.  Global 
differences  such  as  background  illumination,  unequal  intensity  ranges,  and  variable 
contrast  can  impair  the  registration  process  significantly.  To  minimize  the  effect  of 
theses  factors,  the  images  axe  first  quantized  such  that  all  intensity  values  fall  in  the 
range  0  to  255. 

3.3  Global  Alignment 

Although  care  is  taken  during  the  mammography  session  to  ensure  proper 
positioning  of  the  breast,  there  is  very  little  that  can  be  done  to  establish  the  same 
same  alignment  between  two  different  mammograms.  The  registration  system  will 
simply  have  to  be  prepared  to  deal  with  this  problem  in  the  course  of  the  alignment 
procedure. 
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Figure  4  Mammograms  of  the  same  breast  taken  one  year  apart.  Global  alignment 

differences  are  resolved  by  applying  an  affine  transformation  to  the  image 
in  the  upper  left  panel,  (a)  Previous  screening  (b)  Current  screening  (c) 
first  image  after  affine  adjustment  (d)  second  image  unaltered 


The  images  in  Figure  4  are  mammograms  taken  of  the  same  breast  in  two 
different  years.  Although  they  show  the  same  structural  features,  the  breast  on  the 
right  (b)  is  positioned  lower  relative  to  the  breast  on  the  left  and  is  possibly  enlarged. 
It  is  clear  that  the  two  images  must  be  more  closely  aligned  before  attempting  to 
correct  regional  deformations  (e.g.  compression  and  expansion  due  to  differences  in 
patient  positioning  pressure  applied  by  the  imaging  equipment)  . 

3.3.1  Key  Factors.  Analytical  techniques,  such  as  least  squares  fitting, 
offer  one  approach  to  resolving  alignment  differences  (50).  Although  these  methods 
can  perform  quite  well  under  most  circumstances,  they  can  generate  unacceptable 
solutions  if  the  differences  between  the  two  images  are  too  great.  They  are  especially 
sensitive  to  outlier  measurements — those  falling  far  outside  the  acceptable  range  of 
values.  However,  it  is  exactly  these  types  of ’’messy”  situations  on  which  evolutionary 
approaches  thrive.  As  described  in  the  previous  chapter,  genetic  algorithms  excel  at 
searching  through  landscapes  where  rules  are  few  and  far  between  (27,  30). 

The  approach  taken  here  is  to  model  global  (mis)alignment  as  an  affine  trans¬ 
formation  problem.  The  rationale  for  this  decision  is  straightforward:  attempting 
to  resolve  differences  due  to  distortion  cannot  be  done  in  a  meaningful  way  if  the 
images  suffer  from  translational,  rotational,  and  scale  differences.  Once  these  linear 
differences  have  been  resolved,  the  nonlinear  distortions  become  manageable. 

3.3.2  Implementation.  The  global  alignment  subsystem  uses  a  genetic 
algorithm  to  search  for  an  optimal  set  of  transformation  parameters.  The  parameters 
required  to  specify  the  affine  transformation 
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determine  the  transformation  are  sx,sy,9x,9y,  dx  and  dy.  The  four  parameters 
of  the  rotation  matrix,  sx,sy,9x,9y,  specify  the  amount  of  rotation  and  scaling  that 
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will  be  applied,  while  the  dx  and  dy  terms  represent  the  amount  of  translation  along 
the  x  and  y  axes,  respectively,  required. 

A  family  of  25  chromosomes,  each  encoding  the  six  affine  parameters,  are 
allowed  to  evolve  using  crossover  and  mutation  operators.  At  the  end  of  each  gen¬ 
eration,  the  genetic  algorithm  evaluates  the  fitness  of  each  chromosome  using  a 
sum-of-squared  difference  (SSD)  function.  Specifically,  the  parameters  from  each 
chromosome  are  used  to  create  a  new  image,  and  the  calculated  SSD  between  this 
image  and  the  target  image  determine  the  chromosome’s  fitness.  By  minimizing  this 
simple  SSD  function,  the  genetic  algorithm  converges  to  a  specified  error  goal  in 
less  than  20  generations  (in  most  cases).  Table  1  details  how  the  affine  parameters 
are  encoded  in  each  chromosome,  and  Table  2  summarizes  the  genetic  algorithm 
parameters  execution  used  by  the  alignment  module. 


Genetic  Algorithm  Chromosome  String 

Parameter 

Bits 

Description 

dx 

5 

Translation  in  Vertical  Direction 

dy 

5 

Translation  in  Horizontal  Direction 

sx 

4 

Scaling  Ratio  for  Vertical  Direction 

Sy 

4 

Scaling  Ratio  for  Horizontal  Direction 

ex 

4 

Rotation  Angle  X 

dy 

4 

Rotation  Angle  Y 

Table  1  The  six  affine  parameters  encoded  as  genes  in  each  chromosome. 


Genetic  Algorithm  Parameters 

Parameter 

Value 

Chromosome  Length  (bits) 

26 

Population  Size 

25 

Mutation  Rate 

0.25 

Crossover  Rate 

0.9 

Table  2  Genetic  algorithm  execution  parameters. 
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3-4  Multiscale  Decomposition 

This  registration  system  relies  on  a  hierarchy  of  images  to  achieve  the  most  ac¬ 
curate  registration  result  while  minimizing  the  overall  processing  requirement.  Since 
intensity  differences  between  image  segments  can  be  substantial,  more  than  one 
level  will  generally  be  required.  The  general  rule  of  thumb  states  that  the  inter¬ 
frame  velocity  should  not  exceed  one  pixel  per  frame  for  an  accurate  optical  flow 
measurement  (10).  This  system  allows  a  maximum  of  six  levels,  depending  upon 
the  degree  of  variation  between  the  two  images.  Figure  5  depicts  the  basic  pyramid 
structure. 

The  higher  levels  are  all  created  by  smoothing  and  decimation  operations.  To 
obtain  a  new  image  at  level  A:  +  1,  the  image  at  level  k  is  convolved  with  the  2-D 
Gaussian  smoothing  operator  over  a  3  x  3  window  using  the  equation 

2/)  =  4E  H  I(x  +  i,y  +  (13) 

iV  i=-ij=-i 

where 


n  =  '52I(x+hy+o)  (14) 

i=-l  j=- 1 

and  a  is  a  scale  constant.  The  smoothed  image  is  then  decimated  by  a  factor 
of  two  in  each  dimension,  resulting  in  an  image  having  one  fourth  the  area  of  the 
image  one  level  below.  Figures  6  and  7  depict  a  five-tier  pyramid  generated  from 
the  example  source  and  target  images,  respectively.  Each  level  is  magnified  to  full 
base  size  for  ease  of  viewing.  The  axes  show  the  actual  dimensions  of  each  frame. 

The  full  power  of  a  hierarchical  representation  can  be  appreciated  if  we  view  it 
as  the  corresponding  Laplacian  pyramid.  Burt  (10)  developed  the  Laplacian  pyramid 
scheme  to  improve  the  transmission  speed  of  images.  We  can  convert  the  Gaussian 
pyramid  described  above  to  its  Laplacian  counterpart  by  recursively  subtracting  from 
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Multiresolution  Pyramid 


Figure  5  Basic  multiresolution  pyramid  architecture.  Levels  two  and  three  are 
created  by  smoothing  and  decimating  by  a  factor  of  two  the  level  below 
them. 
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the  Gaussian  image  at  level  k  from  the  one  on  the  level  fc  + 1.  Figures  8  and  9  show 
the  Laplacian  pyramids  corresponding  to  the  Gaussian  pyramids  in  Figures  6  and 
7,  respectively. 

Projection  operators  are  needed  to  transform  pixel  values  to  the  appropriate 
values  on  adjacent  levels.  This  system  uses  the  averaging  operator 


2  1 
4  2 
2  1 


and  interpolation  operator 
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1  2 


1 
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1 


(15) 


(16) 


employed  by  Hwang  and  Lee  (33)  in  their  three-tier  Gaussian  pyramid  imple¬ 
mentation. 


3.5  Segmentation 

After  the  two  images  have  been  brought  into  coarse  alignment  by  minimising 
the  translational,  rotational,  and  scale  differences  between  them  and  the  multiresolu¬ 
tion  pyramid  has  been  constructed,  it  useful  to  segment  the  images.  In  this  system, 
two  types  of  segmentation  are  implemented:  breast  structure  extraction  and  internal 
structure  partitioning.  The  first  procedure  extracts  the  region  of  the  image  belonging 
to  the  breast  image  away  from  the  surrounding  background.  This  step  ensures  that 
later  processing  stages  are  matching  breast  tissue  rather  than  attempting  to  match 
pixels  belonging  to  the  background  region. 
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3.5.1  Rationale.  As  discussed  in  Chapter  II,  the  three  most  common 
segmentation  philosophies  are  based  on  either  edge  detection,  region  growing  or 
thresholding  techniques.  The  method  selected  for  this  system  belongs  to  the  second 
group  of  methods  because  it  segments  the  image  by  clustering  similar  pixels  into 
homogeneous  partitions.  Since  the  segmentation  activity  should  operate  without 
the  need  for  outside  intervention,  an  unsupervised  method  would  seem  like  the  best 
choice.  The  Kohonen  self-organizing  feature  map  (SOM)  is  a  neural  network  based  on 
unsupervised  learning  (38).  The  SOM  uses  vector  quantization  to  discover  patterns 
in  the  input  data.  The  Kohonen  SOM  is  unique  among  neural  networks  in  that 
it  creates  a  topographical  map  of  the  input  space.  In  other  words,  patterns  close 
together  in  the  input  space  will  remain  close  to  one  another  in  the  feature  map  as 
well.  The  SOM  map  has  many  of  the  same  characteristics  of  maps  found  in  biological 
information  processing  systems. 

This  thesis  extends  the  previous  method  by  introducing  a  multiresolution  com¬ 
ponent.  Each  image  is  transformed  into  a  deck  of  images  smoothed  by  a  different 
Gaussian  scale  factor.  The  pattern  vector  associated  with  a  given  image  pixel  con¬ 
sists  of  elements  representing  the  pixel  intensity  at  each  scale. 


Kohonen  Extraction  Network 

Parameter 

Value 

Pattern  Vector  Length  (pixels) 

25 

Competitive  Layer  Nodes 

2 

Iterations 

10000 

Update  Mode 

Batch 

Table  3  Kohonen  extraction  network  parameters. 

3.5.2  Breast  Extraction.  Lee  (42)  demonstrated  how  a  Kohonen  self¬ 
organizing  feature  map  (SOM)  can  effectively  partition  a  mammogram  using  a  com¬ 
petitive  learning  paradigm.  First,  a  preprocessing  stage  extracts  small  pixel  neigh¬ 
borhoods  by  sliding  a  5  x  5  window  in  raster  fashion  from  left  to  right,  top  to  bottom. 
These  neighbor  pixel  regions  are  small  enough  to  characterize  homogeneous  areas, 
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Kohonen  Segmentation  Nei 

twork 

Parameter 

Value 

Pattern  Vector  Length  (pixels) 

5 

Competitive  Layer  Nodes 

25 

Iterations 

10000 

Update  Mode 

Batch 

Table  4  Kohonen  segmentation  network  parameters. 

yet  large  enough  to  ensure  the  data  is  statistically  valid.  Figure  10  demonstrates  the 
extraction  capability  of  this  segmentation  technique.  Each  of  the  two  nodes  in  the 
competitive  layer  of  the  network  encode  either  is-breast  or  is-background  based  on 
the  neighbor  intensity  values. 

3.5.3  Breast  Segmentation.  This  thesis  proposes  a  novel  technique  for  seg¬ 
menting  the  internal  structure  of  mammogram  images.  As  with  the  breast  extraction 
procedure  described  above,  this  method  also  uses  a  Kohonen  neural  network.  This 
method,  however,  differs  from  the  other  in  that  it  uses  information  from  the  multires¬ 
olution  pyramid  generated  in  the  previous  processing  step.  Specifically,  an  n-level 
image  deck  is  created  in  which  each  level  is  smoothed  using  a  successively  larger  scale 
parameter  (typically  following  a  geometric  progression).  Then  an  n-dimensional  fea¬ 
ture  vector  composed  of  pixel  values  from  each  level  is  defined  for  each  pixel  in  the 
original  image.  Since  each  successive  Gaussian  image  is  a  low-pass  version  of  the  im¬ 
age  below  it,  pixel  values  at  coarser  levels  will  be  more  resistant  to  noise  than  those 
on  lower  levels.  A  segmentation  based  on  multiple  pixel  values  should  in  principle 
produce  a  partition  that  identifies  more  or  less  homogeneous  regions  in  the  image: 
the  addition  of  more  levels  to  the  image  deck  will  increase  the  smoothness  of  the  par¬ 
titioning.  Figure  11  shows  the  example  source  and  target  images  after  segmentation 
by  the  Kohonen  neural  segmentation  procedure.  Figures  12  through  16  compare  the 
performance  of  the  proposed  Kohonen-based  segmenter  with  the  well-known  fuzzy 
c-means  clustering  algorithm  using  successively  more  pyramid  levels.  Fuzzy  c-means 
clustering  is  described  in  detail  in  (29,  35). 
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Gaussian  Pyramid  [level  two] 


Gaussian  Pyramid  [level  three] 


(c)  (d) 

Figure  6  Gaussian  pyramid  generated  for  source  image  in  Figure  3  [cr  =  1.5].  (a) 
Base  image  decimated  X  4  (b)  second  tier  decimated  X  8  (c)  third  tier 
decimated  X  16  (d)  fourth  tier  decimated  X  32. 


Figure  7  Gaussian  pyramid  generated  for  target  image  Figure  3  [cr  =  1.5].  (a) 
Base  image  decimated  X  4  (b)  second  tier  decimated  X  8  (c)  third  tier 
decimated  X  16  (d)  fourth  tier  decimated  X  32. 
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Lapladan  Pyramid  [level  three] 
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Lapladan  Pyramid  [level  four] 

Lapladan  Pyramid  [level  five] 

Figure  8  Laplacian  pyramid  generated  for  source  image  Figure  3  [a  =  1.5]  •  (a) 
Base  image  decimated  X  4  (b)  second  tier  decimated  X  8  (c)  third  tier 
decimated  X  16  (d)  fourth  tier  decimated  X  32. 


Kohonen  Network  Generated  Binary  Mask 
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(a) 


Kohonen  Network  Segmented  Breast 
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Kohonen  Network  Generated  Binary  Mask 
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(b) 

Kohonen  Network  Segmented  Breast 
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(c)  (d) 

Figure  10  Binary  masks  and  segmented  breast  images  for  two  different  screenings 
[generated  by  Kohonen  Neural  Network],  (a)  Binary  mask  for  previous 
mammogram  (b)  binary  mask  for  current  mammogram  (c)  previous 
breast  structure  extracted  from  background  (d)  current  breast  structure 
extracted  from  background. 
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(c)  (d) 

Figure  11  Mammograms  segmented  using  a  Kohonen  Neural  Network  with  25 
nodes  in  the  competitive  layer.  Segmentation  shows  variations  in  image 
intensity  corresponding  to  regions  of  different  density  (a)  Mammogram 
from  previous  screening  (b)  mammogram  from  current  screening  (c)  seg¬ 
mented  version  of  image  in  upper  left  panel  (d)  segmented  version  of 
image  in  upper  right  panel. 
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(c)  (d) 

Figure  12  Comparison  of  Kohonen  and  Fuzzy  C-Means  (FCM)  segmentation 
methods.  Source  and  target  mammogram  images  are  segmented  into 
8  regions  using  single  level,  (a)  Source  image  segmented  with  Kohonen 
method  (b)  target  image  segmented  with  Kohonen  method  (c)  source 
image  segmented  with  FCM  method  (d)  target  image  segmented  with 
FCM. 
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Figure  13  Comparison  of  Kohonen  and  Fuzzy  C-Means  (FCM)  segmentation 
methods.  Source  and  target  mammogram  images  are  segmented  into 
8  regions  using  two-level  pyramid,  (a)  Source  image  segmented  with 
Kohonen  method  (b)  target  image  segmented  with  Kohonen  method 
(c)  source  image  segmented  with  FCM  method  (d)  target  image  seg¬ 
mented  with  FCM. 
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Figure  14  Comparison  of  Kohonen  and  Fuzzy  C-Means  (FCM)  segmentation 
methods.  Source  and  target  mammogram  images  are  segmented  into 
8  regions  using  three-level  pyramid,  (a)  Source  image  segmented  with 
Kohonen  method  (b)  target  image  segmented  with  Kohonen  method  (c) 
source  image  segmented  with  FCM  method  (d)  target  image  segmented 
with  FCM. 
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Figure  15  Comparison  of  Kohonen  and  Fuzzy  C-Means  (FCM)  segmentation 
methods.  Source  and  target  mammogram  images  are  segmented  into 
8  regions  using  four-level  pyramid,  (a)  Source  image  segmented  with 
Kohonen  method  (b)  target  image  segmented  with  Kohonen  method 
(c)  source  image  segmented  with  FCM  method  (d)  target  image  seg¬ 
mented  with  FCM. 
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Figure  16  Comparison  of  Kohonen  and  Fuzzy  C-Means  (FCM)  segmentation 
methods.  Source  and  target  mammogram  images  are  segmented  into 
8  regions  using  five-level  pyramid,  (a)  Source  image  segmented  with 
Kohonen  method  (b)  target  image  segmented  with  Kohonen  method 
(c)  source  image  segmented  with  FCM  method  (d)  target  image  seg¬ 
mented  with  FCM. 


3.6  Optical  Flow  Computation 

After  the  two  images  have  been  grossly  aligned,  expanded  into  a  multiscale 
representation,  and  segmented  into  homogeneous  regions,  the  optical  flow  between 
the  two  images  can  be  more  reliably  estimated.  The  approach  taken  in  this  thesis 
is  based  on  a  multiscale,  cooperative  computation  paradigm.  Since  we  are  dealing 
with  unknown,  sometimes  significant,  inter-image  differences,  the  multiresolution 
approach  should  be  better  able  to  cope  with  these  challenges  than  its  single-level 
model  counterpart.  The  information  provided  by  the  Laplacian  pyramid  and  that 
obtained  by  residual  errors  further  enhance  the  effectiveness  of  a  hierarchical  strategy. 

3.6.1  Theoretical  Background.  As  described  in  the  previous  chapter,  es¬ 
timating  the  optical  flow  vector  field  can  be  considered  a  constrained  optimization 
problem  having  the  general  form 

$  —  f  f  ( Exu  4-  EyV  +  Et )2  +  ck2( it2  +  Uy  +  v2  +  Vy)dxdy  (17) 

J  J  Image 

where  a  represents  the  Lagrange  multiplier  constant.  The  spatio-temporal 
partial  derivatives  Ex,  Ey,  and  Et  are  approximated  using  the  four-point  differencing 
method  described  in  Horn  and  Schunck’s  original  work  (31).  This  thesis  adopts  the 
solution  method  developed  by  Hwang  and  Lee  (33)  as  outlined  below. 

Since  we  are  using  a  hierarchical  approach,  the  additional  constraint  a2(w2  + 
ul+vl+vl)dxdy  can  be  taken  to  be  the  difference  between  the  velocity  vector  at  the 
current  level  k  and  the  velocity  vector  projected  from  an  adjacent  level  (i.e.  A:  —  1 
or  k  + 1).  Thus,  we  can  model  this  constraint  as  a  minimization  of 

C(uk,  vk)  =  f({(uproj  -  ukf  +  (vproj  -  vk)2))  (18) 

where  (uk,vk)  denote  the  velocity  vectors  at  level  k  and  (upr0J,  upr0J)  are  the 
velocity  vectors  projected  from  an  adjacent  level.  The  minimization  equation  thus 
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becomes 


E(uk,  vk)  =  [  [  (. Exu  +  Evv  +  Etf  +  a2{C{uk,  vk)}2)dxdy  (19) 

J  J Image 

which  in  turn  can  be  formulated  as  a  set  of  Euler-Lagrange  equations  of  the  form 

{a2  +  (Ek}uk  +  EkEkvk  =  a2uproj  -  EkEk  (20) 

EkEkuk  +  {a2  +  (Ek}vk  =  a2vproj-EkEk  (21) 

by  setting  E(uk,  vk)  equal  to  0. 

Hwang  and  Lee  use  the  iterative  Gauss-Seidel  relaxation  technique  to  solve  the 
minimization  equation.  The  iterative  update  equations  are 


«*  =  U0  (22) 

#*  =  (23) 

where 

P  =  Ekuproj  +  Ekvproj  +  Ek  (24) 

D  =  a2  +  (Ek)2  +  (Ek)2  (25) 

The  update  variables  uk  and  vk  are  initialized  to  upr°i  and  vpro\  the  optical 
flow  vectors  computed  at  the  previous  level,  to  provide  a  good  starting  point  for 
further  refinement. 

3. 6.2  Algorithm  Dynamics.  The  coarse- to-fine  computational  strategy  used 
by  this  system  facilitates  the  progressive  refinement  of  an  initial  guess.  In  Hwang 
and  Lee’s  (33)  implementation,  the  optical  flow  calculation  from  the  next  level  above 
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serves  as  the  regularization  parameter  required  to  make  the  problem  well-posed.  This 
design  underlies  the  superior  performance  of  the  multiscale  approach  over  Horn  and 
Schunck’s  (31)  single-scale  approach.  The  proposed  system  can  also  implement  true 
multigrid  processing  by  projecting  the  residual  from  one  level  to  the  level  above  it 
where  it  can  be  used  to  speed  convergence  (19).  Figure  18  depicts  an  example  of  the 
residual  error  between  the  optical  flow  computed  at  two  different  levels  of  the  image 
pyramid.  It  represents  the  high  frequency  component  of  the  computed  flow  field. 
Figure  17  shows  the  final  optical  flow  fields  computed  at  each  level  of  a  four-level 
image  pyramid. 

3. 1  Image  Warping 

The  optical  flow  estimate  can  be  used  directly  to  create  a  warped  version  of 
the  source  image.  This  non-linear  transformation  can  be  accomplished  using  the 
equation 


Imagewarped(i,  j)  —  Imagegonrce(i  Au^j  At?).  (26) 

Since  the  transformed  image  coordinates  (t  —  Au,j  —  Av)  will  be  real-valued 
and  therefore  fall  outside  the  standard  integer  grid,  the  actual  intensity  values  must 
be  interpolated  from  the  source  image.  This  system  uses  a  standard  bilinear  inter¬ 
polation. 

3.8  Evaluation 

The  primary  evaluation  criterion  is  the  similarity  measure  between  the  trans¬ 
formed  image  and  the  target  image.  Although  many  matching  routines  use  some 
type  of  correlation  measure  to  estimate  the  similarity  between  two  images  or  image 
regions,  it  is  not  used  by  this  system  because  it  is  unreliable  when  any  significant 
amount  of  non-rigid  motion  is  present.  Instead,  the  optical  flow  estimate  itself  is 
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used  as  a  measure  of  the  remaining  difference  between  the  two  images.  Since  the 
original  motivation  behind  using  optical  flow  was  to  have  it  estimate  the  inter-image 
disparity,  it  seems  reasonable  that  the  optical  flow  residual  would  likewise  be  a  good 
estimate  of  any  remaining  intensity  difference.  The  lower  the  post-alignment  optical 
flow,  the  lower  the  dissimilarity  (i.e.  better  match)  between  the  two  image  segments. 

Comparison  of  the  two  corresponding  Laplacian  images  offers  a  qualitative 
sense  of  how  well  the  images  match.  As  Barron  et  al.  (3)  point  out  in  their 
evaluation  of  optical  flow  methods,  the  Laplacian  pyramid  ’’helps  to  enhance  image 
structure,  such  as  edges,  that  is  often  thought  to  be  important”.  Whereas  such 
features  are  normally  subjected  to  blurring  in  a  Gaussian  pyramid,  they  can  become 
more  prominent  in  the  Laplacian  pyramid.  If  similar  structural  features  can  be  found 
in  the  warped  and  target  images,  we  can  better  evaluate  the  success  of  the  matching 
procedure. 

3.9  Summary 

This  chapter  presented  a  mammogram  registration  technique  based  on  opti¬ 
cal  flow  estimates  in  a  multiresolution  framework.  After  preprocessing,  the  images 
undergo  coarse  alignment  to  eliminate  major  positional  and  orientation  differences 
between  them.  The  two  breast  structures  are  then  segmented  from  the  background 
image  for  improved  matching  accuracy  and  computational  performance.  Next,  a 
Gaussian  pyramid  of  four  to  five  levels,  depending  on  the  degree  of  disparity  between 
the  images,  is  constructed.  A  multiscale  algorithm  then  computes  the  apparent  dif¬ 
ference  between  the  images  using  an  iterative  coarse-to-fine  approach.  Optical  flow 
estimates  generated  at  a  given  level  become  the  starting  point  for  higher  resolution 
estimates  at  the  next  level  down  the  pyramid.  The  computed  velocity  field  is  used 
to  transform  the  source  image  into  one  more  closely  resembling  the  target  image. 
Finally,  the  disparity  between  the  warped  image  and  the  target  is  evaluated. 
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(c)  (d) 

Figure  17  Optical  Flow  Pyramid:  Decomposition  of  the  optical  flow  estimate  into 
four  levels.  The  calculated  flow  at  level  fc+1  becomes  the  initial  estimate 
for  refining  the  estimate  at  level  k.  (a)  Bottom  level  (finest  resolution) 
(b)  second  level  (c)  third  level  (d)  top  most  level  (coarsest  resolution). 
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IV.  Results 


This  chapter  discusses  results  from  applying  the  proposed  mammogram  registration 
system  to  synthetic  and  real  images. 

4-1  System  Testing 

It’s  important  to  evaluate  the  accuracy  of  the  proposed  registration  method 
before  relying  on  matching  results  using  real  mammograms.  Since  it  is  difficult  to 
precisely  determine  the  expected  results,  the  performance  of  two  different  optical 
flow  estimation  techniques  will  be  analyzed.  Single-level  estimation  will  be  tested 
using  the  Horn  and  Schunck  (31)  method,  while  the  multi-resolution  approach  will 
be  represented  by  Hwang  and  Lee’s  (33)  Gaussian  pyramid  scheme.  Although  not 
exhaustive,  the  tests  are  designed  to  evaluate  how  well  the  proposed  method  can  be 
expected  to  perform  on  real  x-ray  images. 

4-1.1  Image  Preparation.  A  small  mammogram  region  possessing  signif¬ 
icant  texture  variation  was  selected  to  be  the  source  image.  The  image  patch  is 
presmoothed  using  the  Gaussian  filter  described  in  Chapter  III  with  a  4  x  4  kernel 
and  a  standard  deviation  of  1.5  pixels.  Image  intensities  were  then  quantized  to  fall 
in  the  range  0  to  255. 

4-1.2  Test  Case  Generation.  Seven  different  tests  were  conducted  on  image 
pairs  having  known  differences.  The  first  two  test  cases  evaluate  the  ability  of  the 
optical  flow  module  to  recover  translational  motion.  Two  target  images  are  tested: 
one  shifted  down  by  five  pixels  relative  to  the  source  image,  and  another  shifted  to 
the  right  again  by  five  pixels.  For  the  third  test,  the  source  image  is  rotated  about 
the  center  pixel.  Since  the  registration  method  will  be  required  to  perform  nonlinear 
mapping,  the  ability  to  detect  rotational  deformations  is  critical  to  the  success  of 
this  system.  A  similar  test  assesses  how  well  the  system  can  track  movement  of  a 
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edge  contour.  This  test  will  simulate  matching  dissimilar  skin  edges  separated  by 
short  distances.  Finally,  the  reaction  to  finding  an  object  on  the  target  image  not 
present  on  the  source  image  is  investigated.  This  capability  is  important  to  both 
image  mapping  and  subsequent  tumor  detection. 

A  final  test  was  constructed  to  evaluate  how  well  the  system  can  capture  the 
actual  optical  flow  between  two  images.  This  test  was  designed  to  mimic  the  degree 
of  complexity  the  system  would  be  expected  to  handle  in  an  operational  setting. 
The  source  image  in  this  case  was  an  actual  mammogram  downsampled  by  a  factor 
of  16  to  simplify  analysis.  This  source  image  was  then  subjected  to  a  known  affine 
transformation  using  the  technique  discussed  in  the  global  alignment  section.  The 
system  then  had  to  compute  the  flow  between  the  original  image  (source)  and  the 
warped  image  (target).  Ideally,  the  computed  flow  would  be  very  similar  to  the  flow 
defined  under  the  known  affine  transformation. 

4-1.3  Test  Case  Execution.  The  regularization  parameter  a  was  set  to  10, 
and  the  target  error  rate  for  convergence  was  set  to  0.05,  for  all  test  cases.  The  title 
of  each  velocity  flow  diagram  includes  the  total  number  of  iterations  required  for 
convergence. 


Test  Case  Specifications 

Test  Case 

Description 

Parameters 

Horizontal  Translation 

shift  of  target  relative  to  source 

5  pixels 

Vertical  Translation 

shift  of  target  relative  to  source 

5  pixels 

Rotation 

rotation  of  source  about  the  center 

10  degrees 

Translation  and  Rotation 

horizontal  shift  and  rotation 

5  pixels;  10  degrees 

Contour  Tracking 

shift  of  curved  target  boundary 

5  pixels 

Foreign  Object  Insertion 

superposition  of  a  box  onto  source 

10  X  10  pixels 

Table  5  Specifications  for  the  six  validation  test  cases. 


4.1.4  Test  Evaluation.  The  results  of  the  six  tests  are  displayed  in  Fig¬ 
ures  19  through  24.  The  most  apparent  difference 
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In  general,  the  two  methods  produced  remarkably  similar  results.  Two  main 
observations  can  be  made  about  the  comparative  performance  of  the  two  methods. 
First,  Horn  and  Schunck’s  gradient-based  method  invariably  requires  almost  twice  as 
many  iterations  to  converge  as  Hwang  and  Lee’s  multiresolution  approach.  Second, 
the  Horn  and  Schunck  method  is  not  quite  as  consistent  as  Hwang  and  Lee’s  when 
compared  to  the  expected  outcome.  For  example,  almost  all  vector  arrows  in  the 
translation  cases  are  parallel  in  the  Hwang  and  Lee  method,  whereas  a  significant 
percentage  show  distinct  deviation  in  the  Horn  and  Schunck  method.  In  the  contour 
tracking  test,  Horn  and  Schunck’s  method  provided  good  estimates  at  boundary 
regions  only,  while  its  competitor  generated  a  consistent  velocity  field  across  the 
entire  image.  The  combined  accuracy  and  efficiency  favor  a  multiresolution  approach 
for  a  system  as  computationally  intensive  as  the  one  proposed  in  this  thesis. 

The  preferred  multiresolution  method  was  used  in  executing  the  last  test  case 
more  closely  simulating  actual  operating  conditions.  Figure  25  presents  the  results 
of  this  test.  The  computed  optical  flow  is  very  similar  to  the  true  flow.  The  main 
differences  occur  in  regions  of  high  homogeneity.  Gradient-based  methods,  such  as 
the  one  used  here,  generally  find  such  homogeneous  areas  ambiguous  and  difficult  to 
compute  accurately. 

4-2  Mammogram  Registration 

4.2.1  Experimental  Format.  The  primary  focus  of  this  thesis  is  on  the 
registration  of  mammograms  taken  of  the  same  breast  at  different  times.  Therefore, 
registration  experiments  conducted  in  this  section  are  restricted  to  available  temporal 
sequences  of  images.  To  keep  the  analysis  focused  and  easy  to  follow,  this  section 
adopts  a  case  study  approach.  The  two  cases  presented  register  a  time  sequence 
of  four  mammograms  for  the  same  women.  The  first  case  (see  Figure  26)  aligns 
four  craniocaudal  views  for  the  right  breast,  while  the  second  case  (see  Figure  35) 
aligns  four  mediolateral  oblique  views  for  the  same  breast.  The  mammograms  for 
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Figure  19  Comparison  of  two  optical  flow  computation  methods  for  image  transla¬ 

tion  test  case,  (a)  Original  image  (b)  original  image  shifted  to  the  right 
by  5  pixels  (c)  Horn  and  Schunck  gradient  method  (d)  Hwang  and  Lee 
multiresolution  method. 
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Figure  20  Comparison  of  two  optical  flow  computation  methods  for  image  trans¬ 
lation  test  case,  (a)  Original  image  (b)  original  image  shifted  down  by 
5  pixels  (c)  Horn  and  Schunck  gradient  method  (d)  Hwang  and  Lee 
multiresolution  method. 


Figure  21  Comparison  of  two  optical  flow  computation  methods  for  image  rotation 
test  case,  (a)  Original  image  (b)  original  image  rotated  counterclockwise 
by  10  degrees  (c)  Horn  and  Schunck  gradient  method  (d)  Hwang  and 
Lee  multiresolution  method. 
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Figure  22  Comparison  of  two  optical  flow  computation  methods  for  image  transla¬ 
tion  and  rotation  test  case,  (a)  Original  image  (b)  original  image  shifted 
down  by  5  pixels  (c)  Horn  and  Schunck  gradient  method  (d)  Hwang  and 
Lee  multiresolution  method. 


54 


Figure  23  Comparison  of  two  optical  flow  computation  methods  for  image  with 
box  superimposed,  (a)  Original  image  (b)  original  image  shifted  down 
by  5  pixels  (c)  Horn  and  Schunck  gradient  method  (d)  Hwang  and  Lee 
multiresolution  method. 
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Figure  24  Comparison  of  two  optical  flow  computation  methods  for  image  contour 
tracking,  (a)  Original  image  (b)  original  image  shifted  down  and  to  the 
right  by  10  pixels  (c)  Horn  and  Schunck  gradient  method  (d)  Hwang 
and  Lee  multiresolution  method. 
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Optical  Flow  Held  [Theoretical] 
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Optical  Flow  Field  [Computed] 


(c)  (d) 

Figure  25  Comparison  of  computed  optical  flow  with  true  optical  flow,  (a)  Source 
image  (b)  target  image  (c)  actual  optical  flow  field  (d)  computed  optical 
flow  field  using  Hwang  and  Lee  multiresolution  method. 
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each  registration  case  are  arranged  from  (a)  to  (d)  by  date  of  screening  in  ascending 
order.  Figures  27  and  28  are  examples  of  the  Gaussian  and  Laplacian  pyramids 
constructed  for  each  of  the  images  in  preparation  for  the  alignment  process.  All 
images  have  been  enhanced  using  histogram  equalization  to  make  internal  structure 
easier  to  identify. 

4-2.2  Results  Format.  The  results  of  each  registration  set  are  displayed  in 
a  four-panel  graphic.  The  upper  left  and  right  panels  show  the  source  and  target 
images,  respectively.  The  bottom  left  panel  shows  the  final  optical  flow  field  for 
level  three  of  the  Gaussian  pyramid.  This  level  was  chosen  primarily  because  it  is 
easy  to  visualize  on  the  printed  page.  Figure  17  shows  what  all  four  optical  flow 
levels  look  like  for  an  example  registration  session.  The  Laplacian  image  for  the 
first  pyramid  level  is  included  to  help  validate  the  results  of  the  alignment.  Since 
Laplacian  images  typically  enhance  structural  features  at  the  appropriate  scale,  they 
are  useful  in  matching  features  between  two  images  and  thus  evaluating  registration 
performance. 
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(c)  (d) 

Figure  26  Temporal  sequence  of  four  mammograms  taken  during  different  screen¬ 
ings.  These  images  are  used  for  individual  two-image  matches  in  the 
next  three  case  studies.  All  images  are  craniocaudal  views  of  the  left 
breast.  The  sequence  (a)-(d)  represents  the  order  in  which  the  mammo¬ 
grams  were  taken,  with  (a)  being  the  earliest  and  (d)  being  the  latest. 


4-2.3  Registration  Case  1. 


4-2. 3.1  Key  Factors.  Figure  26  displays  the  four  mammograms 
comprising  registration  case  1.  All  images  are  craniocaudal  (CC)  views.  The  global 
alignment  stage  is  not  required  for  these  four  images  due  to  the  absence  of  any 
significant  translation  or  rotation  difference  between  them. 

4-2.3. 2  Registration  Results.  The  registration  system  has  no  diffi¬ 
culty  lining  up  the  breast  outline  in  any  of  the  three  image  sets.  Also,  the  internal 
structures  (based  on  intensity  characteristics  only)  correspond  well  between  images. 
Case  1.1  shows  the  most  difference  in  the  upper  third  of  the  images.  The  computed 
optical  flow  field  is  particularly  prominent  in  that  region  and  weak  everywhere  else. 
The  flow  field  in  case  1.2  clearly  demonstrates  the  ability  of  the  system  to  detect 
translational  motion  between  two  images.  The  smooth  vertical  flow  is  consistent  with 
a  downward  shift  in  the  target  image  as  can  be  readily  by  comparing  the  two  breast 
outlines.  Case  1.2  presents  a  challenge  because  it  is  not  at  all  clear  whether  the  dark 
patch  in  the  lower  left  quadrant  is  due  to  tissue  compression  of  an  artifact  of  poor 
imaging.  Since  the  registration  procedure  cannot  differentiate  either,  it  attempts  to 
map  the  dark  region  in  the  source  and  target  images.  The  resulting  warped  image 
displayed  in  panel  (d)  matches  up  well  in  the  upper  half  of  the  images.  In  each 
case,  the  system  computed  flow  fields  consisting  of  both  linear  (translational)  and 
non-linear  (distortion)  components.  It  is  precisely  this  ability  to  deal  with  regional 
distortion  differences,  such  as  compression  and  expansion,  that  makes  the  approach 
so  appealing  and  potentially  useful. 
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Figure  29  Case  1.1  Registration  of  two  mammograms  taken  at  different  times 
[craniocaudal  view].  The  optical  flow  field  shown  corresponds  to  level 
2  of  the  Gaussian  pyramid  to  facilitate  viewing  (a)  Mammogram  from 
previous  screening  (b)  mammogram  from  current  screening  (c)  optical 
flow  field  calculated  from  source  image  [in  upper  left  panel]  to  target 
image  [upper  right  panel]  (d)  mammogram  from  current  screening. 


a006b00bjnam  [Lapladan  Pyramid]  a006b!0u.mam  [Lapladan  Pyramid] 
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Figure  30  Case  1.1  Lapladan  images  of  the  warped  and  target  images,  (a)  Source 
image  (b)  target  image  (a)  warped  image  from  source  (b)  target  image. 
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(c) 


(d) 


Figure  32 


Case  1.2  Laplacian  images  of  the  warped  and  target  images,  (a)  Source 
image  (b)  target  image  (a)  warped  image  from  source  (b)  target  image. 
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4-2-4  Registration  Case  2. 


4-2-4-1  Key  Factors.  Figure  35  portrays  the  four  mammograms  used 
for  registration  case  2.  Whereas  the  previous  sets  involved  the  craniocaudal  view, 
case  2  deals  with  images  taken  from  the  medial-lateral  oblique  (MLO)  view.  In 
most  MLO  views,  the  pectoral  wall  can  be  clearly  seen  and  can  serve  as  an  excellent 
registration  landmark.  A  significant  vertical  shift  is  noticeable  between  the  source 
and  target  images.  This  translation  was  not  adjusted  by  the  alignment  module  in 
order  to  assess  the  ability  of  the  system  to  include  such  a  shift  in  the  computed 
optical  flow  field. 

4-2. 4-2  Registration  Results.  In  case  2.1,  the  pectoral  wall  protrudes 
deeper  in  the  target  image  than  it  does  in  the  source  image.  The  computed  flow  field 
easily  detects  this  strong  intensity  shift  as  seen  in  the  strong  translational  response 
in  the  upper  left  of  the  flow  field.  The  obvious  global  misalignment  in  case  2.1 
was  purposely  not  handled  by  the  alignment  module  to  demonstrate  the  ability  of 
the  multiresolution  algorithm  to  resolve  large  intensity  shifts  between  images.  The 
computed  flow  displays  a  strong  vertical  component  in  an  attempt  to  compensate  for 
the  significant  downward  translation  from  source  to  target  image.  Closer  inspection 
of  the  flow  field  shows  that  the  algorithm  is  able  to  deal  with  both  the  relative 
vertical  shift  and  inward  compression  as  seen  in  the  center  right  portion  of  the  flow 
field. 

Case  2.2  is  instructive  because  it  shows  how  the  system  reacts  when  a  signifi¬ 
cant  portion  of  the  source-target  pair  exhibit  only  minor  differences,  while  elsewhere 
large  differences  are  visible.  As  expected,  the  computed  flow  field  is  weak  between 
the  highly  similar  regions  and  strong  between  the  regions  exhibiting  large  intensity 
differences.  Case  2.3  requires  a  large  vertical  downward  shift  and  displacement  of 
the  pectoral  wall  to  the  left. 
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Figure  33  Case  1.3  Registration  of  two  mammograms  taken  at  different  times 
[craniocaudal  view].  The  optical  flow  field  shown  corresponds  to  level 
2  of  the  Gaussian  pyramid  to  facilitate  viewing  (a)  Mammogram  from 
previous  screening  (b)  mammogram  from  current  screening  (c)  optical 
flow  field  calculated  from  source  image  [in  upper  left  panel]  to  target 
image  [upper  right  panel]  (d)  mammogram  from  current  screening. 
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Figure  35  Temporal  sequence  of  four  mammograms  taken  during  different  screen¬ 
ings.  These  images  are  used  for  individual  two-image  matches  in  the 
next  three  case  studies.  All  images  are  medial-lateral  oblique  views  of 
the  left  breast.  The  sequence  (a)-(d)  represents  the  order  in  which  the 
mammograms  were  taken,  with  (a)  being  the  earliest  and  (d)  being  the 
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Figure  36  Case  2.1  Registration  of  two  mammograms  taken  at  different  times 
[mediolateral  oblique  view].  The  optical  flow  field  shown  corresponds 
to  level  2  of  the  Gaussian  pyramid  to  facilitate  viewing  (a)  Mammo¬ 
gram  from  previous  screening  (b)  mammogram  from  current  screening 
(c)  optical  flow  field  calculated  from  source  image  [in  upper  left  panel]  to 
target  image  [upper  right  panel]  (d)  mammogram  from  current  screen- 
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Figure  38  Case  2.2  Registration  of  two  mammograms  taken  at  different  times 
[mediolateral  oblique  view].  The  optical  flow  field  shown  corresponds 
to  level  2  of  the  Gaussian  pyramid  to  facilitate  viewing  (a)  Mammo¬ 
gram  from  previous  screening  (b)  mammogram  from  current  screening 
(c)  optical  flow  field  calculated  from  source  image  [in  upper  left  panel]  to 
target  image  [upper  right  panel]  (d)  mammogram  from  current  screen- 


Figure  39  Case  2.2  Laplacian  images  of  the  warped  and  target  images,  (a)  Source 
image  (b)  target  image  (a)  warped  image  from  source  (b)  target  image. 


Table  6  summarizes  results  for  both  case  1  and  case  2  registration  sets.  The 
similarity  metric  between  images  is  a  sum  of  squared  difference  taken  between  each 
pixel  intensity  value  in  the  first  image  and  the  corresponding  pixel  intensity  value  in 
the  second  image.  Difference  1  measures  the  total  intensity  difference  between  the 
source  and  target  image  pair.  The  second  measures  the  difference  between  the  source 
and  the  warped  image  obtained  from  the  source  image,  while  the  third  measures  the 
difference  between  the  warped  and  target  images.  The  pattern  is  clear:  in  each 
case,  the  difference  between  warped  and  target  image  is  less  than  half  the  difference 
between  source  and  target  image  pairs.  This  result  shows  that  the  match  between 
warped  and  target  images  is  better  than  between  the  original  source-target  image 
pair. 

4-3  Summary 

This  chapter  presents  results  obtained  from  test  case  simulations  and  the  regis¬ 
tration  of  historical  mammogram  sets.  The  optical  flow  produced  during  the  testing 
phase  corresponded  very  well  with  expected  test  results.  The  multiresolution  ap¬ 
proach  in  general  proved  superior  to  the  single  image  method  by  achieving  a  higher 
degree  of  consistency  and  executing  much  faster.  Acceptable  optical  flow  estimates 
were  obtained  after  one  relaxation  sweep  (approximately  50  iterations). 

Results  obtained  from  registering  actual  mammogram  images  were  consistent 
with  expected  results. 
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Figure  40  Case  2.3  Registration  of  two  mammograms  taken  at  different  times 
[mediolateral  oblique  view].  The  optical  flow  field  shown  corresponds 
to  level  2  of  the  Gaussian  pyramid  to  facilitate  viewing  (a)  Mammo¬ 
gram  from  previous  screening  (b)  mammogram  from  current  screening 
(c)  optical  flow  field  calculated  from  source  image  [in  upper  left  panel]  to 
target  image  [upper  right  panel]  (d)  mammogram  from  current  screen¬ 
ing. 
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Figure  41  Case  2.3  Laplacian  images  of  the  warped  and  target  images,  (a)  Source 
image  (b)  target  image  (a)  warped  image  from  source  (b)  target  image. 
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Registration  Results  Summary 

Case/Set 

Iterations 

Difference  1 

Difference  2 

Difference  3 

1.1 

55 

6.04 

4.04 

2.9 

1.2 

56 

4.85 

3.63 

2.01 

1.3 

51 

5.25 

3.81 

2.17 

2.1 

38 

5.00 

3.04 

2.30 

2.2 

34 

3.31 

2.15 

1.40 

2.3 

50 

3.54 

2.70 

1.33 

Table  6  Summary  of  registration  statistics.  The  first  the  difference  between  origi¬ 
nal  image  pair.  The  second  is  the  between  source  image  and  warped  image. 
The  third  is  the  difference  between  warped  image  and  target  image.  All 
difference  measures  are  sum  of  squares. 
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V.  Conclusion 


This  chapter  summarizes  the  proposed  registration  method,  discusses  contributions, 
and  suggests  possible  avenues  for  future  research. 

5.1  Summary 

Breast  cancer  is  lethal  disease,  indiscriminately  afflicting  nearly  two  hundred 
thousand  women  in  this  country  each  year.  Early  detection  of  suspected  cancerous 
growth  is  an  important  weapon  in  the  fight  against  this  deadly  killer.  Today’s  mam¬ 
mography  screening  procedures  provide  an  effective  means  toward  this  end  but  in 
general  fail  to  achieve  the  desired  level  of  detection  accuracy.  Radiologists  typically 
compare  mammograms  taken  at  different  times  to  improve  the  likelihood  that  a  can¬ 
cer  will  be  noticed.  Although  computer-aided  detection  (CADx)  are  being  developed 
to  help  the  radiologist  make  a  more  informed  decision,  current  systems  are  generally 
limited  to  single-image  analysis. 

The  image  registration  system  developed  in  this  thesis  attempts  to  bridge  the 
gap  between  the  needs  of  the  radiologist  during  a  mammography  screening  and  the 
capabilities  of  an  automated  assistant.  Its  primary  purpose  is  to  align  two  mam¬ 
mograms  of  the  same  breast  so  that  the  radiologist  can  make  a  direct  comparison 
between  a  time  sequence  of  images.  Since  it  is  often  very  difficult  to  establish  struc¬ 
tural  features  in  the  image  pair  to  facilitate  a  mapping  between  them,  this  system 
is  designed  not  to  rely  on  the  availability  of  such  landmarks.  Instead,  the  proposed 
system  compares  regional  intensity  values  and  calculates  the  apparent  motion  be¬ 
tween  source  and  target  image.  This  computed  optical  flow  field  then  forms  the  basis 
for  a  non-linear  mapping  between  the  two  images.  This  transformation  results  in  a 
warped  version  of  the  original  source  image  which  can  be  directly  compared  with 
the  target  image.  Finally,  the  underlying  architecture  of  the  registration  system  is 
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based  on  a  multiresolution  computational  model,  significantly  enhancing  both  the 
accuracy  and  efficiency  of  the  matching  process. 

The  proposed  system  incorporates  features  to  provide  a  robust  registration 
result.  Since  many  image  pairs  exhibit  significant  differences  in  the  positioning 
and  orientation  of  the  breast  structure,  a  global  alignment  module  is  included  to 
correct  for  such  differences  prior  to  computing  the  mapping  between  the  images. 
This  module  relies  on  a  genetic  algorithm  to  quickly  and  accurately  search  for  an 
optimal  global  alignment.  It  requires  no  user  intervention  to  perform  the  alignment 
task.  A  second  module  is  provided  to  perform  critical  segmentation  duties.  This 
segmentation  subsystem  uses  a  neural  network  to  partition  each  image  into  regions 
of  high  homogeneity.  The  module  is  designed  to  carry  out  two  different  segmentation 
tasks:  1)  extraction  of  the  breast  structure  from  the  surrounding  background,  and 
2)  partitioning  on  the  breast  structure  itself  into  regions  of  similar  intensity.  The 
segmented  image  can  often  help  the  radiologist  interpret  an  image  and  facilitate 
identification  of  different  tissue  types. 

The  initial  performance  of  the  system  was  evaluated  on  two  sets  of  images 
taken  of  the  same  breast.  Each  set  of  images  represented  a  historical  sequence  of 
mammograms  taken  at  four  different  screenings.  One  set  consisted  of  four  cran- 
iocaudal  views,  while  the  other  presented  the  complementary  mediolateral  oblique 
views.  Several  test  cases  where  constructed  and  executed  to  evaluate  system  per¬ 
formance.  Then  a  mapping  between  pairs  of  images  (source  and  target)  in  each  set 
was  computed  to  obtain  a  warped  source  image.  The  system  performed  well  even  in 
the  face  of  major  differences  between  images.  In  all  cases  the  breast  contour  of  the 
warped  and  target  images  matched  up  very  closely.  The  system  demonstrated  its 
proficiency  at  coping  with  distortion  such  as  compression  and  expansion  by  generat¬ 
ing  a  non-linear  optical  flow  field.  In  all  cases,  the  sum-squared  intensity  difference 
between  warped  and  target  image  pairs  was  always  considerably  less  than  the  original 
difference  between  source  and  target  images. 
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5.2  Contributions 


As  a  part  of  the  effort  to  develop  the  mammogram  registration  system,  three 
specific  contributions  can  be  highlighted: 

1.  Application  of  optical  flow  and  multiresolution  computation  to  solve  the  mam¬ 
mogram  correspondence  problem 

2.  Application  of  a  genetic  algorithm  to  find  an  optimal  set  of  parameters  for  an 
affine  transformation 

3.  Development  of  a  segmentation  method  based  on  the  Kohonen  neural  network 
and  multiresolution  image  decomposition 

5.3  Recommendations 

The  research  conducted  in  preparing  this  thesis  does  not  purport  to  solve  the 
mammogram  correspondence  problem.  The  task  of  bringing  two  images  into  com¬ 
mon  alignment — especially  images  having  significant  differences — is  a  complex  and 
perhaps  unsolvable  problem.  The  work  presented  here  attempts  only  to  provide  some 
fresh  ideas  that  might  stimulate  even  better  approaches.  The  system  developed  in 
this  thesis  should  be  improved  and  expanded  to  directly  address  cancer  detection. 
Further  insight  into  the  problem  might  lead  to  a  system  capable  of  flagging  poten¬ 
tial  tissue  abnormalities  as  its  computes  the  optical  flow  field  between  two  images. 
Much  creative  work  remains  toward  the  development  of  an  intelligent,  comprehensive 
approach  to  image  comparison  and  abnormality  detection. 
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13.  ABSTRACT  (Maximum  200  words ) 

Breast  cancer  is  second  only  to  lung  cancer  as  the  most  prevalent  form  of  cancer  to  afflict  women — remaining 
the  leading  cause  of  cancer  death  in  women  between  the  ages  of  40  and  55.  Mammography  is  a  potent  weapon  in 
the  fight  against  this  lethal  disease,  due  in  large  part  to  its  widespread  availability  and  low  cost.  Despite  the  fact 
that  mammography  can  detect  small  lesions  as  early  as  two  years  before  they  become  palpable  on  physical  exam, 
between  10  and  30  percent  of  cancerous  lesions  go  undetected  during  evaluation  by  the  radiologist.  One  approach 
to  improving  detection  rates  involves  comparing  mammograms  of  the  same  breast  from  successive  years.  Since 
most  forms  of  breast  cancer  develop  slowly,  multiple-view  techniques  might  be  able  to  detect  subtle  changes 
indicative  of  cancerous  growth.  This  thesis  proposes  a  computer-aided  system  designed  to  bring  two  images 
into  correspondence,  or  alignment,  so  that  they  can  be  compared  and  evaluated  for  possible  abnormalities.  The 
system  estimates  a  mapping  between  two  images  by  calculating  the  optical  flow,  or  apparent  intensity  change, 
between  a  source  and  target  mammogram.  The  efficiency  of  the  proposed  registration  system  is  enhanced  by 
utilizing  a  multiresolution  approach  whereby  images  are  compared  at  more  than  one  scale.  In  contrast  to  other 
registration  attempts  which  match  sets  of  morphological  features,  this  system  does  not  require  the  identification 
of  any  control  points  at  all.  This  advantage  permits  the  system  to  perform  well  even  when  the  two  images  differ 
significantly  from  one  another.  Preliminary  results  suggest  the  potential  usefulness  of  this  system  as  part  of  a 
clinical  computer-aided  detection  (CADx)  system. 
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